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We propose a model for increasing liquid saturation in a granular packing which can account for 
liquid redistribution at saturation levels beyond the well-studied capillary bridge regime. The model 
is capable of resolving and combining capillary bridges, menisci and fully saturated pores to form 
local liquid clusters of any shape. They can exchange volume due to the local Laplace pressure 
gradient via a liquid film on the surfaces of grains. Local instabilities like Haines jumps trigger 
the discontinuous evolution of the liquid front. The applicability of the model is demonstrated and 
compared to benchmark experiments on the level of individual liquid structures as well as on larger 
systems. 

PACS numbers: 47.15.gm, 47.55.-t, 47.56.+r, 68.03.-g, 68.15.+e 


I. INTRODUCTION 

The study of liquid distribution and transport in unsat¬ 
urated porous media, such as granular packings, has long 
been a topic of intense cross-disciplinary research. The 
mechanical properties of granular materials are signifi¬ 
cantly affected by their liquid saturation m as liquid 
clusters emerge and grow in size. While the liquid bridge 
regime at low saturation is well-studied numerically and 
experimentally Em, the extension of models to higher 
liquid contents remains a demanding task. A better un¬ 
derstanding of this regime is however crucial for solving 
a number of open problems in science and engineering 
including e.g. rainfall-induced slope failures, oil recovery 
or flow in porous media, just to name a few. 

Early works in this field employed the so-called ideal 
soil model, consisting of uniform solid spheres in a regu¬ 
lar packing 00 . In the 1950s, network models emerged, 
originally proposed by Fatt [9]. They represent only the 
pore space by sites (pore bodies) of arbitrary shape and 
position, interconnected by bonds (pore throats) to form 
a network PEHE]. Network topologies are either ob¬ 
tained from thin section analysis mi, tomographic mea¬ 
surements [15] , from fundamental statistical assumptions 
[131 or directly from packed particle configurations [T6] . 
Even though the geometrical representation of the pore 
space is significantly simplified, important information 
can be obtained for two-phase flows such as relative per¬ 
meability m, trapped immiscible fluid [13], drainage 
and imbibition m, or capillary pressure - saturation re¬ 
lations [18]. A major limitation of pore network models 
lies in the static representation of the pore space with re¬ 
spect to deformations, where network parameters would 
constantly change. 

To cope with such challenges, a more geometrically de¬ 
tailed representation of the pore space and of the fluid 
interfaces is needed. While the pore network is still de¬ 
rived from a random sphere packing that could evolve in 
time, details on liquid interfaces like liquid bridges and 
menisci of constant curvature |T are considered [1T9H2T] . 
In these grain-based fluid invasion models of three di¬ 


mensional porous media, the fluid front is driven by local 
instabilities and correctly reproduces drainage and imbi¬ 
bition experiments. 

With the availability of advanced microtomography, 
a rich variety of liquid clusters was found [22] [23] . The 
number and the size of observed liquid clusters was shown 
to strongly depend on the saturation level. Inspired by 
these findings, we propose a model that explicitly consid¬ 
ers all possible liquid morphologies on the pore scale in 
a sphere packing. Our aim is to be able to simulate the 
entire range of saturation levels from the dry state, via 
cluster formation and growth to the fully saturated state. 
Even though we use a discrete element model (DEM) 
with spherical particles to simulate the random packing, 
we fix all particle positions for the fluid simulation. Hence 
we describe the first step on the way to a general model 
that combines two-phase flow and deformation. 

First, in chapter [H] we give a comprehensive description 
of the model with respect to implemented liquid struc¬ 
tures and their composition into larger clusters along 
with geometrical stability criteria for their growth and 
decay. In the following chapter m we provide a brief 
overview on the simulation procedure. We show appli¬ 
cations of our model on two validation experiments; one 
on the pore scale for the trimer formation and decay and 
the other one for the evolution of cluster distributions 
in large systems (chapter |IV| ) . Finally we show a typical 
simulation of the cluster evolution for the case of liquid 
injection at a singular point, before we summarize main 
results and draw conclusions (chapter [V]). In Appendix [B] 
we describe in detail the important pressure update al¬ 
gorithm for the volume controlled simulations. 


II. COMPONENTS OF THE GRAIN SCALE 
MODEL 

Our aim is to represent both the pore space and the 
liquid structures as close to their real shape as possible, 
accepting drawbacks in numerical performance. Hence 
the pore-throat network is constructed based on the exact 
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geometrical positions of a previously simulated granular 
packing. Liquid clusters are represented as a combina¬ 
tion of elementary units such as liquid bridges, menisci 
or entirely filled pore bodies. Not only units can form 
higher geometrical configurations such as trimers, pen- 
tamers, tetrahedral clusters and higher ones, but they 
can also decompose into elementary ones. This evolu¬ 
tion is determined by stability criteria discussed in this 
chapter as well. 



FIG. 1. (Color online) Single tetrahedral cell from a trian¬ 
gulated particle packing on sphere centers Pi, i — 1,.., 4 with 
radius R. S 23 denotes the contact distance of particle 2 with 
respect to particle 3 and 023 the respective opening angle. 

We extract the pore space from the initial particle 
packing given by the positions of the sphere centers 
and their radius R via Delaunay triangulation following 
Ref. [16]. Hence the entire sample volume is subdivided 
into single tetrahedra like the one shown in Fig. [l] The 
void space in each tetrahedron is called the pore body 
while the cutting areas of the pore body with the respec¬ 
tive faces of the tetrahedron form the four pore throats 
of a cell. This way, a pore network is extracted directly 
from the topology of the sphere packing. Pore bodies can 
be empty, partially filled with liquid separated by menisci 
or entirely filled by liquid. In the discussed model we as¬ 
sume that the pressure inside the liquid phase is smaller 
than the gas pressure. 


A. Representation of Liquid Structures 

When pore bodies are filled with liquid, they are con¬ 
sidered as one of the three building units of the liquid 
clusters, while liquid bridges and menisci shown in Fig. [2] 
are the other two. A liquid bridge is located between 
two grains (Fig. |5Ja)). The pressure difference between 
gas and liquid phase A P = Pu qu id — Pgas is due to sur¬ 
face tension 7 and well described by the Young-Laplace 


equation [24] 

iP = 7C - 7 (i- + i-), (1) 

where the surface curvature C = 1/Ri + I/R 2 is de¬ 
fined by the principal radii of curvature R\ and R 2 . This 
equation can be solved numerically for a liquid bridge 
with boundary conditions defined by the contact angle 
0. Since solving the Young-Laplace equation is a nu¬ 
merically expensive and difficult task, we interpolate the 
capillary pressure P and the volume V of a liquid bridge 
from tabulated values. These were calculated by Sem- 
prebon et al. [25| using the numerical energy minimiza¬ 
tion method of the software Surface Evolver [26]. Both 
pressure P and volume V are functions of the separation 
distance S between the grains, the contact angle 6 for 
solid-liquid interfaces and the filling angle (3 of the liquid 
bridge. Note that 0 is kept constant at 5° throughout 
this work. 

A key parameter for a liquid bridge and its pressure is 
the separation distance Sij between the grains i and j. 
Note that this distance remains constant during the fluid 
simulation, however capillary bridges can rupture due to 
liquid outflow. The dimensionless rupture distance S c 
of the liquid bridge (in units of the particle radius R) 
is related to the dimensionless volume V and the con¬ 
tact angle 0 through the empirical expression derived by 
Willett et al. m- 

s c ~ (1 + ±e)(W + vVvio). (2) 

An important assumption of our model, motivated by 
ideas of Haines 0, is that the liquid-air interface between 
three grains, called meniscus is of spherical shape with 
constant curvature (Fig. [ 2 ]). A cross section perpendic¬ 
ular to the pore throat is shown in Fig. [2jc) . This ap¬ 
proximation is in satisfying agreement with experimen¬ 
tal observations [22] [23] . The centers of the four possi¬ 
ble menisci inside a tetrahedral cell are located on the 
normal of each pore throat through the circumcenter of 
the respective face. The exact position on this normal 
can be calculated once the contact angle with the grains 
0 and the meniscus radius R men are known, following 
the method proposed by Gladkikh [ 20 ] (see Fig. [ 2 |. The 
pressure drop between gas and liquid is again determined 
by the curvature of the interface 1 /Rmen via the Young- 
Laplace equation (Eq. [TJ. Due to the spherical shape of 
the meniscus the negative pressure drop is given by 

A-P = 7 • 7 T“ • (3) 

amen 

In reality single menisci cannot exist, but always occur in 
combination with associated liquid bridges forming one 
liquid body with equal Laplace pressure and common 
surface (see Fig. [ 2 | [ 22 , [23j [28]. We adopt the assump¬ 
tion of equal Laplace pressure within the liquid body, 
determined by the curvature of the meniscus. Note that 
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FIG. 2. (Color online) (a) Inter particle liquid bridge with separation distance S 12 and filling angle shown in the cross section 
through the centers of the particles Pi — P 2 . Ri ,2 denote the principal radii of the curvature, (b) Trimer: Meniscus in the pore 
throat which is build by three grains with centers Pi,2,3- The concave meniscus (green) has one contact point (61,2,3) with 
each neighboring particle. Additionally, connected liquid bridges (blue) are shown, (c) Enlarged cut through P 1 -C 1 -K showing 
the meniscus position in the pore throat. The point K is located between 62 and 63 on the meniscus, see (b). The point O 
is the center of the meniscus. Z614 and W 23 denote the liquid bridges between the particles Pi and Pa respectively P 2 and P 3 . 
The hatched area shows the cut through the cylinder which approximates the meniscus volume. Figure (c) is based on [ 20 | . 


this pressure can fall below the values for stable bridges, 
resulting in a liquid body with less than three associ¬ 
ated liquid bridges. In small clusters menisci with less 
than three associated liquid bridges have not been ob¬ 
served [22 . Therefore, we allow for menisci with a re¬ 
duced number of bridges when clusters have more than 
one filled pore body. This assumption is required to en¬ 
able clusters to fill also those regions of the pore space 
where inter particle separations are large. 



FIG. 3 . (Color online) Small liquid clusters: (a) Pentamer 
and (b) tetrahedral cluster. 

A main innovation of our approach is the possibility 
for the formation of local liquid clusters which can 
evolve inside the granular material. They are composed 
of the three basic units introduced above: filled pore 
body, liquid bridge, and meniscus. In principle every 
liquid body with the exception of a single liquid bridge 


can be considered as a cluster. The smallest possible 
cluster is called trimer (Fig. [2|)) [23]. In our model it 
is bounded by two menisci that are located in neighbor¬ 
ing tetrahedra on both sides of the common pore throat 
and their shared liquid bridges. Higher order clusters 
like pentamers are formed when two trimers share one 
liquid bridge (Fig.[3^). The smallest cluster with a filled 
pore body is a tetrahedral cluster with the liquid body 
bounded by four menisci and six associated liquid bridges 

(Fig#)- 

The control variable in our model is the volume, thus 
its correct calculation for the liquid clusters is essential 
for any simulation. Remember that inside a liquid clus¬ 
ter, pressure is homogeneous and given by the radius of 
its menisci Rmen • In general, the volume of a cluster V c 
with Nimiy imbibed pore bodies and N men menisci can be 
written as 

N irnb AT men 

Tc(-^men) = ^ ^ Vpore,i T ^ ^ I^nen,j (-^men) 5 (4) 

i=0 j=0 

where V pore ^ denotes the volume of the imbibed pore 
body i and Vm en j the volume of the meniscus j includ¬ 
ing the associated liquid bridges of this meniscus. Note 
that liquid bridges in a cluster are always associated with 
one of its menisci creating a common liquid body. The 
volume of the filled pore is calculated by subtracting the 
partial volumes of the four particles contained within the 
tetrahedral cell from the volume of this cell. 

The volume beneath a meniscus is approximated by 
the volume of a cylindrical body V cy i coaxial to the nor- 
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mal of the pore throat, see Fig. p[c). Its radius R cy i 
equals the circumradius R c irc = P\T minus the particle 
radius R, and its upper surface is bounded by the menis¬ 
cus. Note that if the meniscus bounds a filled pore body, 
the meniscus can intersect the pore throat. In this case 
the volume of the pore body beneath the throat plane 
must be reduced by the volume bounded between the 
meniscus and the throat pore plane. As mentioned above, 
to calculate V men ^ we also take into account volumes of 
the liquid bridges V- b associated with the meniscus i: 

m 

imen(-^men) k cyl (R 

men ) + (0.5+e) j2 v I b ( R men )• (5) 

i=0 

The geometrical correction parameter e accounts for the 
volume excess of a real meniscus with the connected liq¬ 
uid bridges, compared to our approximation in which V cy i 
underestimates the volume bounded by the meniscus (see 
Fig-i c)). Depending on the number of the connected 
liquid bridges the index m can take values between 0 and 
3. If a meniscus shares a liquid bridge with an adjacent 
meniscus (e.g. for a pentamer, Fig. [3^i) the number of 
connected liquid bridges for one of the menisci reduces 
by one. With this rule we assure that the bridge volumes 
are not counted twice. Note that for every meniscus only 
half of each liquid bridge volume is considered, since the 
other half is located in the opposite triangulation cell if 
its pore body is not saturated (see the liquid bridge lb \4 
between the particles Pi and P 4 in Fig. ie). 

B. Evolution of Liquid Structures 

Since we assume all particles to be fixed in space at 
this stage, liquid structures can evolve when their volume 
changes due to condensation or evaporation at gas-liquid 
interfaces or the injection of liquid at distinct positions. 
Inside of a liquid body the flux is instantaneous, while 
transport between bodies i and j, sharing the same par¬ 
ticle occurs through a liquid film on the particle surface. 
The relevance of liquid transport through this wetting 
layer was experimentally proven by Lukyanov et al. m ■ 
Under the assumption of stationary flow in the film, the 
volume flux is considered to be proportional to the local 
capillary pressure gradient A P = Pi — Pj. The dimen¬ 
sionless flux into a liquid structure i is then given by the 
sum of all volume fluxes between i and the connected 
structures j : 

TD Ni 

= (6) 

7 3 =0 

where Ni is the number of all liquid structures connected 
to i. The dimensionless conductance coefficient fiij must 
include effects associated with the geometrical distance 
between structures, the number of structures connected 
to one particle etc [25, 30i. However, for simplicity we 
assume this coefficient to be fiij = 0.01 for all calculations 


in this paper. The value of fi^ defines the time scale 
of the pressure equilibration between connected liquid 
structures. Note that this liquid transport mechanism is 
in particular important at low saturation levels. 

The propagation of the liquid interface in the mate¬ 
rial is triggered by local instabilities. Hence we model 
interface evolution as a discontinuous process with in¬ 
stantaneous jumps from one stable configuration to the 
next one depending on the local volume, if the corre¬ 
sponding pore body is drained or imbibed. An initially 
stable interface can become unstable due to curvature 
changes, when the volume of the liquid body increases or 
decreases. We implemented five geometrical instability 
criteria similar to the ones described by Motealleh et al. 
[ 2 Tj . The first four result in growth while the last crite¬ 
rion accounts for drainage resulting in shrinkage of the 
liquid cluster: 

Criterion cl: Pore throats are filled due to coales¬ 
cence of liquid bridges if they touch each other, creating 
a new trimer. The criterion is based on the filling angles 
Pi and fa of the bridges and the opening angle of the 
corresponding pore throat 2a^. If 0.5(/?i + P 2 ) > 2a^ 
the pore throat is imbibed (see Fig. j2jb)). 

Criterion c2: Imbibition of a pore body due to an 
interface instability is triggered when a meniscus touches 
a single liquid bridge inside the tetrahedral cell that was 
up to now not part of the liquid body (Melrose criterion, 
m)- The transition from trimer to tetrahedral cluster 
(Fig. [3^i —^b) is such an example. Using the filling angle 
of the meniscus the configurations become unstable 
when + P/2 > p with the face-edge angle p shown 
in Fig. [21(c) . The calculation is described in detail in 
Ref. [20^ 

Criterion c3: Imbibition of a pore body due to an 
interface instability can also be triggered by the merging 
of two menisci in the same pore body. The criterion 
is a generalization of the Haines imbibition criterion to 
non-zero contact angles. If two menisci centers reach the 
same point, the menisci touch and build a single spherical 
interface. Then the menisci become unstable, disappear 
and the pore body is filled [ 21 . Note that in situations 
where all four possible menisci of a tetrahedral cell exist, 
a gas bubble gets trapped. 

Criterion c4: When a meniscus touches the opposite 
particle, the pore body gets entirely filled [21]. This event 
occurs if the curvature of the meniscus l/R men is small 
enough to allow a touch of the meniscus and the opposite 
grain (the exact value depends on the geometry of the 
particular triangulation cell). Note that this criterion is 
an extension of the one proposed in Ref. [32] from two to 
three dimensions. 

Criterion c5: Decreasing volumes result in increasing 
absolute values of Laplace pressure that can trigger two 
different kinds of drainage: (c5-l) the decay of trimer 
units of clusters and (c5-2) the drainage of pore bod¬ 
ies. Hence the first one describes breakup of a liquid 
body inside a single pore throat. Due to the assumption 
of a meniscus with constant curvature, an intuitive cri- 
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terion of touching menisci on opposite sides of the pore 
throat overestimates the stability of the trimers by about 
20% with respect to experimental observations [ 22 ] . We 
propose a criterion based on a minimum thickness for 
stability: 


HZin > *R, (7) 

with the adaptable drainage parameter k to calculate the 
critical height H™™ of menisci for stability. Best agree¬ 
ment with experiments is obtained for a value of k = 0.15. 
Depending on the opening angle 0^23 (see Fig. [ 6 ]), two or 
three liquid bridges remain. The drainage criterion (c5- 
2) is related to drainage of a single pore body. When the 
center of the meniscus of a neighboring cell touches the 
respective pore throat plane of an entirely saturated cell, 
its pore body is drained. Instantaneously the liquid inter¬ 
face jumps to a new stable position, creating three new 
menisci in the pore throats of the drained pore body, re¬ 
moving the responsible meniscus in the neighboring cell. 
This criterion follows the ideas in Ref. m and is an 
extension of the original one proposed by Haines [7 to 
arbitrary contact angles. Now that the model description 
is complete, we address its simulation. 


III. NUMERICAL IMPLEMENTATION 


Prior to the fluid simulation, we construct the pore 
space, defined by a dense particle packing, using a 
discrete element method (DEM) described in detail in 
Ref. [30] with contact dynamics. We perform a random 
sequential adsorption [331 of equally sized spheres inside 
the sample volume V sarnp i e with periodic boundary con¬ 
ditions. To achieve a desired packing density p, parti¬ 
cle radii R are increased while particles are allowed to 
rearrange. A small amount of liquid is assigned to ev¬ 
ery particle to stay in the pendular regime where only 
liquid bridges exist. Whenever particles contact, a liq¬ 
uid bridge is created; when particles separate beyond S c , 
bridges rupture Q. In the following, the liquid content 
is defined using the total volume of the liquid structures 
Vliquid as W c = Vliquid/Vsampie • Finally, particles are 
fixed at their positions for the remainder of the simula¬ 
tion, liquid bridges are allowed to equilibrate pressure, 
and the sample is subdivided into tetrahedral cells by a 
periodic Delaunay triangulation of particle centers using 
the CGAL package [34] . 

In our implementation the volume is the control vari¬ 
able. Liquid is either injected (removed) at distinct point 
in space or condensed (evaporated) at gas-liquid inter¬ 
faces. This results in an update of all liquid structures to 
make the pressure correspond to the new volume. The 
cluster pressure update algorithm returns the pressure in 
a single liquid cluster after its volume or configuration 
has changed. This algorithm is explained in detail in 
Appendix [B] After the liquid c onten t has changed, all in¬ 
stability criteria (cl-c5) (Sec. IIB) are checked. Before 


pore bodies can be filled (c2-c4), menisci must form, cor¬ 
responding to criterion cl. Hence, first a list of bridges 
fulfilling cl is generated. Then new trimer units are 
formed by these bridges, eliminating all involved bridges 
from the list for further evaluation. Note that whenever 
a trimer is added to an already existing cluster due to 
cl, the entire cluster needs to be updated. The same 
holds if two clusters are merged by a trimer unit. When 
no more bridges fulfill criterion cl, we proceed with the 
other criteria in a similar fashion. This sequential proce¬ 
dure is chosen, since after fulfilling the higher numbered 
criterion, the preceding ones become even more unlikely 
and the liquid structure more stable. For criteria c2- 
c4 also a list of unstable menisci is created, from which 
the most unstable configuration is chosen and processed 
(pressure update in the corresponding cluster). Then, 
the list of unstable menisci is updated. Note that alter¬ 
native sequences of c2-c4 were tested with similar out¬ 
come. When liquid is drained from the system, mainly 
criterion c5 is relevant, however cl-c4 are checked after 
the pressure update for consistency. 

Due to the instabilities, we obtain pressure differences 
in liquid structures that drive liquid transport through 
liquid films (see Sec. IIB, Eq.[ 6 |. After the transport, the 
cluster pressure is recalculated to account for changes in 
volumes of liquid bodies. After the final pressure update, 
the time is incremented by At. In Appendix we show 
the program sequence (Fig. 0 and list all simulation 
parameters (Tab. [I]). 


IV. CHANGING LIQUID CONTENT ON THE 
PORE SCALE 


Trimers are one of the building units for all larger clus¬ 
ters in a particle packing and the most common morphol¬ 
ogy after the liquid bridge for a wide range of saturation 
values [ 22 ] • As they are well studied and rather simple, 
trimers are ideal validation and calibration clusters for 
the numerical model e.g. for trimer creation (criterion 
cl) and decay of single trimers (criterion c5-l). For the 
case of condensation, we study the distribution of cluster 
morphologies as a function of liquid content and compare 
to experiments by Scheel et al. [ 22 ]. We demonstrate the 
applicability of our method for the case of liquid injection 
at singular points until full saturation. 


A. Trimer formation and decay 

A trimer in our model consists of two menisci that are 
located on both sides of a single pore throat and three 
connected liquid bridges (Fig Dp), analogous to observa¬ 
tions with microtomography ^3] . We calculate trimer 
creation by merging of two or three liquid bridges (crite¬ 
rion cl) and their decay (c5-l), and compare to experi¬ 
mental data [22 (see Fig. [4|. 
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FIG. 4. Dimensionless Laplace pressure of a trimer and liq¬ 
uid bridges as a function of their volumes normalized by the 
number of contacts N cont and the grain radius R 3 . Trimers 
form at Pmax and decay at Pmin for a contact angle of 6 = 5°. 
The experimental data denoted by squares with error bars is 
taken from Ref. m ■ The lines show the simulation results. 


opening angles result in a high minimal pressure Pmim 
e.g. for 0^23 > 32.5° no trimers can exist at a Laplace 
pressure Ptrimer < —4.4. For imbibition this implies that 
trimers with particles not being in contact are generated 
at higher overall Laplace pressures compared to contact¬ 
ing ones. 



V/N cont R 3 


To reproduce the experiment we first simulate conden¬ 
sation into liquid structures. For this, a small amount of 
liquid AV con d = 0 . 002 i ? 3 is added to each one of the three 
liquid bridges between the contacting grains in every time 
step of the simulation (see Fig. 11 in the Appendix). The 
liquid bridge pressure P is updated to correspond to the 
changing volume. At a critical dimensionless value Pmax , 
the bridges coalesce and a new trimer is created (Fig. [4|. 
At the same time the liquid pressure drops, since liquid 
is needed to fill the pore throat between the grains, thus 
decreasing the radius of the interface curvature. Then 
the process is inverted with liquid evaporating from the 
surface of the trimer. In this case, a small amount of liq¬ 
uid AV evap = 0.01 R 3 is removed from the trimer in ev¬ 
ery time step while its pressure is adjusted. Eventually, 
the trimer decays into three liquid bridges at the critical 
minimal dimensionless pressure Pmin • We confront our 
model prediction for this process with the experimental 
data in Fig. [4j To obtain good agreement the geometri¬ 
cal correction factor e is set to 0.07 (see Eq. [5| and the 
drainage parameter k for the instability criterion c5-l 

(Eq -0 was set to k =0.15. 

The critical minimal pressure P m in of a trimer strongly 
depends on the separation distance Sij between particles, 
or on the corresponding opening angles aij for a monodis- 
perse particle packing. Since three liquid bridges are re¬ 
quired for a trimer, we assume that a trimer can no longer 
exist if the separation distance Sij is larger than the rup¬ 
ture distance S c for one of the associated liquid bridges 
(see Eq.[2|. In Fig. [5] we show the pressure-volume curves 
for a trimer in a particle configuration in which one of 
the opening angles otij = 0^23 (see Fig. [2jb)) is increased. 
Note that 0^23 = 30° corresponds to the previously stud¬ 
ied case of a closed contact (Fig. [4]). As expected, large 


FIG. 5. Dimensionless Laplace pressure as a function of the 
normalized volume of a trimer for different opening angles 023 
(see Fig. |T]) . Also indicated are the minimal Laplace pressures 
Pmin for which the trimer is still stable. 



<*23 [ ] 

FIG. 6. Trimer decay: critical minimal Laplace pressure as a 
function of the opening angle 023 for trimers. Experimental 
data by Scheel m ■ Three types of trimer decay are observed 
in the experiment: (1) decay into three liquid bridges, (2) 
decay into two liquid bridges, (3) slow transition into two 
liquid bridges. Only configurations with one gap opposite to 
the opening angle are shown. 

We demonstrate the validity of the above assump¬ 
tion by comparing the calculated dimensionless critical 
Laplace pressures P m in as a function of the opening an¬ 
gle c ^23 for single trimers with measurements by Scheel 
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[22] . In the experiment, trimers can decay into three 
(30° < 0^23 < 32°) or two liquid bridges with fast 
(32° < (^23 < 34°) and slow ( 0^23 > 34°) transition, 
as shown in Fig. [ 6 ] In the simulation liquid is evapo¬ 
rated from trimers in a dense packing through reduction 
of their volumes by AV evap in every time step. In Fig. [ 6 ] 
the critical pressure P m ^ n of the single trimer decays is 
shown for configurations with one gap opposite to the 
opening angle 0 ^ 23 . Only decays of trimers are recorded 
in which the other two gaps are less than 10 - 2 P. Our 
results are well within the error bars of the experimen¬ 
tal measurements. We observe two different regimes for 
opening angles 0^3 > 31.7° and 0^3 < 31.7°. The former 
is the result of the rupture of the liquid bridge, while the 
latter is due to the drainage criterion iJ™*™ — with k 
selected as k = 0.15 (see Eq.[7|. The particular value of k 
was chosen to deliver the best agreement with experimen¬ 
tal observations in Figs. [4] and [ 6 ] The critical pressure 
curve was also calculated analytically by considering a 
single trimer in the three particle configuration with a 
varying opening angle 0 ^ 23 . This curve is also shown in 
Fig.© Small deviations of the simulated single trimer 
decays from the analytical solution are due to the finite 
size of AV evap in the simulation. 


B. Morphogenetics of liquid clusters 

Studies on single liquid clusters are important for the 
model validation, however liquid in a random packing 
leads to a multitude of clusters of different topology 
and size. Their statistic was evaluated from microtomo- 
graphic images for a random close packing with density 
p = 0.57 =b 0.01 at different liquid contents W c [ 22 ]. In 
the simulation we increase W c by condensing liquid into 
clusters and single liquid bridges starting with a pendu¬ 
lar regime. A small amount of liquid proportional to the 
surface is added to every structure in each time step. The 
obtained distributions are shown and compared with ex¬ 
perimental ones in Fig. [7] A snapshot of the sample with 
different cluster morphologies at W c = 0.03 is shown in 
Fig- m At low liquid contents only liquid bridges ex¬ 
ist, while larger clusters start to emerge later. Above 
the pendular state, the number of clusters rapidly in¬ 
creases and reaches a maximum which depends on the 
corresponding cluster size. After the maximum of the 
respective morphology is reached, the clusters predomi¬ 
nantly merge and build larger structures. Note that fi¬ 
nally only one single percolating cluster remains. Our 
simulation reproduces the general trend observed in the 
experiment. However, some deviations are present (in 
particular for tetrahedral clusters). One of the main rea¬ 
sons for these deviations, apart from the approximate 
nature of instability criteria cl-c5, is a finite sample size 
of the simulation with only up to 11 single tetrahedral 
clusters being observed (see Fig. | 8 |. 


C. Haines Jumps in Single Cluster Growth 

When liquid is injected into the particle packing at a 
singular point, pressure drops in the single liquid body 
of the emerging cluster occur. Pressure drops originate 
from rapid imbibitions of single pore bodies, so called 
Haines jumps, and have been observed in experiments 
(e.g. Ref. [35]). Pressure drops are due to volume con¬ 
servation with respect to the new configuration, resulting 
in increasing curvature of menisci. Starting from a liq¬ 
uid content in the pendular regime ( W c « 1 %), we inject 
liquid with a constant flux into the center of the sam¬ 
ple. Initially, liquid bridges at the injection point grow 
and eventually coalesce to a trimer that grows until the 
first instability occurs. Then the cell pore body is im¬ 
bibed, creating a tetrahedral cluster. In the following, 
an increasing number of pore bodies is imbibed, result¬ 
ing in a drop of Laplace pressure inside the cluster (see 
Fig# In the beginning, the dimensionless cluster pres¬ 
sure increases and reaches the value of —0.79, followed 
by a sharp drop. Each drop is associated with the local 
instability of filling a new pore body or throat. In Fig. [9] 
we observe decreasing magnitudes of pressure drops with 
increasing cluster size. This can be explained by the fact 
that the amount of liquid needed to fill the next pore 
body in relation to the total cluster volume decreases 
with increasing cluster volume. Note that the pressure 
converges to a approximately constant value. This value 
is given by the minimal possible pressure of a ” critical” 
liquid bridge of the cluster which has a maximal sepa¬ 
ration distance between the particles (see Fig. |5j. The 
morphology of the growing cluster at different time steps 
is shown in Fig. [9] Prior to the first pressure drop, the 
liquid cluster is a trimer (t = 0.15), then a tetrahedral 
cluster (t = 0.25), and later on it grows in all direc¬ 
tions in a compact manner, even with single trimer units 
being connected to it (t = 5.01). The compact shape 
of the cluster is due to the pressure equilibration time 
scale which is set by the conductance coefficient p. For 
p = 0.01 the liquid transport through surface films is 
slow compared to the cluster growth that is triggered by 
the liquid injection rate. Therefore the cluster grows in a 
compact manner mainly through imbibition of pore bod¬ 
ies. The pressure-volume relationship for a trimer like 
the one at the beginning of cluster growth (Fig. 91) has 
been compared with experimental results in Fig. 4] For 
larger clusters no experimental data is available. How¬ 
ever, the convergence to a constant pressure level during 
the cluster growth can be linked to experimental observa¬ 
tions which state that the pressure in equilibrated liquid 
morphologies remains constant in a wide range of satu¬ 
ration levels [23] . 

V. SUMMARY AND CONCLUSIONS 

We propose a model for fluid saturation in random 
packings that can cope with arbitrary liquid contents 





FIG. 7. (Color online) Number of liquid morphologies N scaled with the maximal number of bridges Nbridge in a packing 
with the density p exp = 0.57 d= 0.01 (experiment) respectively p S im = 0.57 (simulation) as a function of the liquid content W c . 
Experimental data from Ref. [22] . 



is able to reproduce all kinds of liquid clusters on the 
grain scale observed in experiments like bridges, trimers, 
pentamers and so on. This was shown on a number of ex¬ 
perimental benchmark problems. We demonstrated that 
by using simple geometrical approximations, fairly good 
agreement with experiments can be obtained at only a 
small portion of effort required to model the true shape of 
the interfaces by for example energy minimization meth¬ 
ods [26] . Since we use volume as a control variable, we 
are capable of calculating pressure drops in single clus¬ 
ters due to Haines jumps, as shown in the example of 
Sec. |IV C| The approach is not limited to static pore 
spaces, and can be used with minor modifications for 
the study of unsaturated deformable granular packings. 
Studies on the interaction of the pore space with the 
pressure field are of particular interest for questions of 
soil behavior. 


FIG. 8. (Color online) Snapshot of the sample with different 
cluster morphologies at W c = 0.03. Only liquid bodies of the 
clusters are shown. Single liquid bridges are omitted for clar¬ 
ity. Green: trimers, magenta: pentamers, grey: heptamers, 
blue: tetrahedral clusters. Large clusters are colored accord¬ 
ing to their volume with a color map ranging from yellow to 
red. 


ranging from dry to full saturation with good accuracy. 
Volume is used as a control variable and pressure is cal¬ 
culated based on the volume conservation for every liquid 
cluster contributing to the liquid content in the porous 
sample. This innovation is of particular importance since 
the volume respectively the water content is much easier 
to access experimentally than the pressure. The model 
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FIG. 9. (Color online) Dimensionless Laplace pressure of a liquid cluster as a function of time, when liquid is injected at a 
constant flow rate at a fixed point with evolving liquid cluster in terms of its liquid body (insets). 


TABLE I. Summary of simulation parameters. 


Label 

Value 

Definition 

N p 

2000 

number of particles 

P 

0.57 

packing density 

R 

1.31 

particle radius 

At 

0.001 - 0.01 

time step 

0 

5° 

contact angle 

7 

1 

surface tension 

P 

0.01 

conductance coefficient 

e 

0.07 

geometrical correction parameter 

K, 

0.15 

drainage parameter for meniscus 


Appendix A: Simulation parameters 



FIG. 10. Calculation of the final meniscus radius Rmen which 
lies between the radius with the minimal volume Vmin and 
the one with a larger volume Vmax here shown for a trimer. 


Simulation parameters are summarized in Tab. [T| 


Appendix B: Pressure update algorithm 


The main aim of the pressure update algorithm is 
to find the new cluster radius Rmen such that the 
Vc(Rmen) ~ Void = 0 is satisfied, where V Q \d denotes the 
cluster volume before and V c after the pressure update. 
This root-finding problem is numerically solved using the 
false position method (Regula Falsi) see Fig. 10 The tar¬ 
get radius of the recalculated cluster menisci radius Rmen 
which corresponds to the volume V c is located between 
Rmin and R ma x • V m in and Vmax denote the cluster vol¬ 
umes calculated with the radius Rmin respectively Rmax • 
The value R m in is the minimal radius for which the re¬ 
calculated cluster is still stable. This lower threshold is 
given by the pressure of the most unstable meniscus tak¬ 
ing into account the associated liquid bridges. An im¬ 
plication of choosing a lower threshold R m in is that in 
some cases the available volume Vm is smaller than the 


minimal volume of the new (recalculated) cluster Vmin'- 
Void < Vmin(Rmin )• This is for example the case when a 
new pore body should be imbibed by liquid from a small 
cluster which has not enough liquid to fill the pore body 
and create a stable cluster interface. Since in this case the 
imbibition event cannot occur, the corresponding insta¬ 
bility is saved in the list of impossible instabilities which 
are not tested again in the current time step. The maxi¬ 
mal radius Rmax can be chosen depending on the cluster 
volume before the volume change or instability event. It 
should be noted here that there is a saturation volume 
for every cluster which corresponds to the case where all 
menisci of the cluster have infinite radius. In real simula¬ 
tions this volume is never reached because an instability 
will appear before. Therefore it is always possible to find 
a proper value for R m ax such that V c < V ma x • Tor de¬ 
creasing cluster volumes, the cluster radius before the 
pressure update is used as R ma x • If the cluster volume 
increases, Rmax is set to a multiple of the cluster radius 
before the pressure update. In the end, the pressure up¬ 
date algorithm delivers the radius of the menisci in the 
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Creation of a particle packing with density p and 
initial water content W c . 

Triangulation via CGAL package. 

for t 

= 0 to t = T end with increment At 


Add/remove liquid to/from structure(s). 

Pressure update in clusters. 

Pressure update in single liquid bridges. 

for i= 0 to i= Nsingle ib. with increment 1 


Overlap of single l.b .? 

Create new trimer 

Check cl for menisci bridges. 

Check criteria c2 -c4. 

Water transport in liquid films. 


Pressure update in clusters: 


for i 

— 0 to i = N c i U s te r with increment 1 


Calculate new radius R of cluster i (see App. B) 

Rmen exists? . V 


Drainage c5 : Sns°taSe memscus 

Drain pore body/throat and calculate 
the new radius R men 

Recalculate all menisci of cluster i with R men 


Pressure update in liquid bridges: 


for i 

= 0 to i = N sing i e i b with increment 1 


Assign new pressure from data base 


Assign new filling angle 


Water transport: 


for k 

= 0 to k = Nu qu id with increment 1 


for i = 0 to i = N str ucture with increment 1 

Volume change with Eq. 6. 

Pressure update in clusters. 

Pressure update in single liquid bridges. 


Checking of criterion cX (X=l,...4): 
unstable = true 


for 

i = 0 to / = N ce n with increment 1 



configurations for cX2^"^F 


Add i to the list of unstable config. 



\ Unstable configurations exist? 


Choose the most unstable configuration 

unstable 

-false 

Imbibe pore/throat and calculat cR m en,new 

Recalculate all menisci of cluster / 
with the radius Rmen,new 



cluster Rmen for which the volume is conserved. 


FIG. 11. Scheme of the simulation progress. 
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